##################################################################
# /* Author: Gautam Nair and Nicholas Sambanis */
# /* Violence Exposure and Ethnic Identification: Evidence from Kashmir */
# /* Appendix Table 12: Refusal Rates by Treatment Condition  */

##################################################################

# change working directory here
# setwd("")

##################################################################
# loading packages
##################################################################
library(Hmisc)
library(foreign)
library(sandwich)
library(lmtest)
library(numDeriv)
library(stargazer)
library(ggplot2)
library(plyr)
library(gridExtra)
library(ri)
library(dplyr)
library(scales)

##################################################################

rm(list=ls())

library(stargazer)

data.working <- read.csv("nei_kashmir_replication_dataset.csv")

dep.var <- c("r_identity_rank_indian_rev",
"r_identity_choice",
"r_kashmir_status",
"r_ind_pak_prefer",
"r_bonus_donation",
"r_ic_index_id_no_ind_pak",
"r_ic_index_id_ind_pak",
"r_protest_peaceful_endorse",
"r_protest_violent_endorse"
)

data.working$group_rand1_rand2_2 <- as.numeric((data.working$group_rand1_rand2==2))
data.working$group_rand1_rand2_3 <- (data.working$group_rand1_rand2==3)
data.working$group_rand1_rand2_4 <- (data.working$group_rand1_rand2==4)
data.working$group_rand1_rand2_5 <- (data.working$group_rand1_rand2==5)
data.working$group_rand1_rand2_6 <- (data.working$group_rand1_rand2==6)
data.working$group_rand1_rand2_7 <- (data.working$group_rand1_rand2==7)
data.working$group_rand1_rand2_8 <- (data.working$group_rand1_rand2==8)

ind.var <- c("group_rand1_rand2_2", "group_rand1_rand2_3" , "group_rand1_rand2_4", "group_rand1_rand2_5", "group_rand1_rand2_6", "group_rand1_rand2_7", "group_rand1_rand2_8")

reg.models <-vector("list", length(dep.var)) 
se.models <-vector("list", length(dep.var))

for(i in 1:length(dep.var)){
	temp.data <- data.working[,c(dep.var[i],ind.var)]
	temp.data <- na.omit(temp.data)
	dep.var[i]
	my.formula <-paste(dep.var[i],'~', paste(ind.var, collapse= ' + '))
	temp.model <- lm(my.formula, data=temp.data)
	temp.se <- coeftest(temp.model, vcov=vcovHC(temp.model,type="HC2"))[,2]
	se.models[[i]] <- temp.se
	reg.models[[i]] <- temp.model
}
reg.models

varlabels <-       c("Violence Exposure", "Economic Growth", "Economic Growth + Violence", "Integrative Institutions", "Integrative Institutions + Violence", "Geographic Priming", "Geographic Priming + Violence", "Constant (Control Group Mean)")
spectitle <- c("Refusal Rates by Treatment Condition")
outputfile <- c("tf_t_12_refusal_rates")
temptype <- c("latex", "text")
tempext <- c(".tex", ".txt")

for(q in 1:length(temptype)){
	temp.output <- paste(outputfile, tempext[q], sep="")

 stargazer(reg.models,se=se.models,
         title= spectitle,
          out=temp.output,
          no.space=TRUE, 
          model.numbers= TRUE,
         notes.append=FALSE,
          align=FALSE,
         covariate.labels=varlabels,
          header= FALSE,
          font.size="small",
          model.names=FALSE,
          digits=3,
         star.char=c("**", "***"),
         star.cutoffs=c(0.05, 0.01),
          omit.stat = c("rsq", "f","adj.rsq"),
          column.labels  = NULL,
          dep.var.caption  = "",
          dep.var.labels.include=FALSE,
         df=FALSE,
         type=temptype[q],
         omit.table.layout="n"
) 
}